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Abstract 

A TVD scheme has been developed and incorporated into an existing time- 
accurate high-resolution Navier-Stokes code. The accuracy and the robustness of 
the resulting solution procedure have been assessed by performing many 
calculations in four different areas: shock tube flows, regular shock reflection, 
supersonic boundary layer, and shock boundary layer interactions. These 
numerical results compare well with corresponding exact solutions or experimental 
data. 


1 


1. Introduction 


Recently, an iterative-implicit diagonally dominant factorization algorithm, 
together with a high order finite-difference scheme, for solving the multi- 
dimensional compressible unsteady Navier-Stokes equations has been developed. 
The important features of this solution algorithm, the finite difference scheme and 
some validated results were reported in Ref. [1], The present work is a continuing 
effort at developing its shock-capturing capability through the use of flux limiters. 
A brief description of the resulting TVD (total variation diminishing) scheme is 
given in Section 2. 

Standard test cases have been carried out to assess the overall accuracy of 
the current code, which embodies the solution algorithm/scheme presented in Ref. 
[1] and the shock-capturing capability developed under the current effort. In 
Section 3, four different shock-tube flows are calculated to demonstrate its 
capability of accurately tracking transient flow discontinuities including shock 
waves and contact discontinuities. In Section 4, the results of an oblique shock 
wave reflection are shown and compared with exact solutions. This demonstrates 
its accuracy in computing shocks that are not aligned with grid lines. To evaluate 
the possible effects of flux limiters on the numerical solution of flows containing 
boundary layers, calculations of a supersonic laminar boundary layer flow have 
been performed, the results are presented in Section 5. Using this viscous flow over 
a flat plate as the starting condition, an oblique shock wave is then introduced to 
examine the interaction of an oblique shock wave with a laminar boundary layer. 
Both the transient development of the interaction and the comparison of the steady- 
state solution between the current calculation and the corresponding experiment 
are discussed in Section 6. 
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2. The Solution Algorithm and The TVD Scheme 


A relatively detailed discussion of the overall solution algorithm and a high 
resolution finite-difference scheme for the inviscid fluxes, termed as the FCTD 
scheme, can be found in Ref. [11. In short, the basis of the solution algorithm is a 
diagonally dominant approximate factorization procedure. The factorization error 
and the timewise linearization error associated with this baseline procedure are 
reduced by performing Newton-type inner iterations at each time step. The 
robustness of the overall algorithm is enhanced by carrying out the temporal 
iterations in pairs to enforce the operational symmetry of the factorization 
procedure. The temporal accuracy is increased to second-order by using three-point 
backward time differencing. The viscous fluxes are evaluated by using the half- 
spacing second-order central differencing scheme. The inviscid fluxes are evaluated 
by the so called FCTD scheme, which is an amended fourth-order central 
differencing scheme with its injected numerical dissipation having the same form as 
the entire dissipative part of the truncation error intrinsic to the third-order-biased 
upwind scheme. Under the current effort, a TVD form of the FCTD scheme has 
been developed to capture flow discontinuities that often occur in many practical 
problems. A brief description of this development follows. 

The two dimensional Navier-Stokes equations in generalized coordinates (%,r|) 
can be written as 


8Q 

dt 


+ ±{e-e u 



= 0 


(1) 


/v 

where Q = Q/ J ; Q = (p,pu,pv,e) T , and J = ^ x r| y -^ y T| x is the metric Jacobian. Here, T 
is the time; p is the fluid density; e is the total internal energy per unit volume; u 
and v are the velocity components in the x and y directions of a Cartesian 
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coordinates system. The transformed inviscid fluxes are denoted by E and F and the 
transformed viscous fluxes are denoted by E v and F v . The specific forms of these 
transformed fluxes are well known and will not be repeated here. Using an 
iterative implicit technique and three point backward time differencing, equation (1) 
leads to 


where 


and 


-Q_5Q + -?J 8£ - + AJ& . 5F V 
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= rhs m 
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/ = 1 , 2 , 3 , ...m 


In the above equations, f denotes an arbitrary quantity, l represents an iteration 
index, and m is an intermediate iteration level between the n-th and (n+l)-th time 
levels. 

The construction of operators approximating the left hand side of Eq. (2) and 
the evaluation of the viscous fluxes in Eq. (3) are the same as those described in Ref. 
[1], The development of a TVD scheme for evaluating the inviscid fluxes in Eq. (3) 

/s 

is based on the FCTD scheme described in Ref. [1], Taking — in Eq. (3) as an 

/s 

example, and dropping ' a * from the flux E, the working formula for FCTD scheme is 
given by 

[FCTD](E)i = (l-p)[FC](E)i + p[TU](E)i (4) 
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where [FC](E)i and [TU](E)i denote the fourth-order central differencing and the 
third-order-biased upwind differencing of at the spatial point (i j) respectively. 

It is noted here that, for convenience, the index j associated with the redirection has 
been dropped in Eq. (4). 

The numerical dissipation of the FCTD finite-difference scheme is essentially 
an infinite series with its elements being the fourth and higher even derivatives of 
the absolute fluxes. The relative amount of added dissipation can be controlled 
through an adjustable parameter p. When P=0, FCTD becomes the fourth-order 
central differencing scheme without numerical dissipation and when p=l, it 
becomes the third-order-biased upwind scheme. 

Following the steps detailed in Ref. [2], a TVD form of the FCTD scheme can 
be constructed, and it has the following form : 


where 


and 
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where E is the total flux, E+ and E" are the split fluxes, A+ is the forward 
differencing operator, and $ is a limiter. In Eqs. (8) and (9), the symbol < > denotes 

the inner product. The expression of Ea is obtained by replacing i with i-1 in the 

above equations. 

Now comes the choice of the limiter and the flux splitting scheme. In the 
present work, Roe’s "superbee" [3] and flux-difference splitting [4] are chosen. Such 
a preference is based on our experience with several Euler and Navier-Stokes 
solutions. Nevertheless, it should be pointed out that no limiter has been found to 
be universally satisfactory, and a drawback of Roe's flux-difference splitting is that 
it may not spread the expansion wave correctly. 

3. Unsteady Shock Tube Problems 

Numerical results for four different shock tube cases have been obtained by 
solving the 2-D Euler equations and compared with the exact 1-D solutions. Table 1 
gives the normalized initial conditions of these standard test cases. For Sod's and 
Lax's problems, the computational domain is covered by 101x21 uniformly 
distributed grid points (Ax=Ay=0.1). For the strong shock and large temperature 
ratio problems, 201x41 points (Ax=Ay=0.05) are used. At the upstream as well as 
the downstream boundaries, flow quantities are consistently over specified by 
invoking the exact solutions. Along the top and bottom boundaries, symmetry 
conditions are imposed. The one-dimensionality of the calculated results has been 
confirmed, and the results along the horizontal centerline of the computational 
domain will be presented. 

The Sod and Lax problems involve only moderate strength shock. For the 
Sod problem, constant At=0.03 is used to advance the solution over a period of 60 
time steps. Results for values of (3 between 0 and 1 have been obtained and 
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compared well with the exact solution. As an example, the solutions obtained with 
[1=1.0 are shown in Fig. 1. For the Lax problem, constant At=0.017 is used to 
advance the solution over a period of 85 time steps. Again, various [1 values have 
been used and typical results for (3=1.0 are presented in Fig. 2. For the large 
temperature ratio problem, fine grids and constant At=0.0051 are used. In this 
case, satisfactory results can only be obtained with [3 > 0.75 . The results shown in 
Fig. 3 are obtained with (3=1.0 over a period of 100 time steps. There are some 
oscillations in the neighborhood of shock and contact discontinuity. 

As mentioned before, an apparent drawback of Roe's flux-difference splitting 
is that it may not spread the expansion wave correctly and thus lead to a non- 
physical 'expansion shock' appearing in the computed flow. Fig. 4 illustrates such a 
case. For this shock tube flow, the pressure ratio is 30, and the density ratio is 24. 
The grid spacings are Ax=Ay=0.05, and a constant At=0.01 is used to advance the 
solution over a period of 90 time steps. The comparison between the exact solution 
and the computed results obtained from the current TVD scheme with (3=1 (i.e., the 
third-order-biased upwind scheme) indicates the appearance of an expansion shock’ 
at x=0. It is noted here that the same 'expansion shock' also occurs for other values 
of [3. When the pressure ratio is increased to the present strong shock level of 500, 
the solution procedure quickly diverges as a consequence of a much stronger, non- 
physical, expansion shock. One way to fix this problem is to add dissipation 
proportional to the strength of the expansion (see e.g. [3]). Under the current effort, 
the fourth-difference constant-coefficient artificial dissipation model [5] is adopted 
for this purpose. Fig. 5 illustrates the solutions obtained with (3=1.0 and o=1.5, 
where a is the value of the dissipation coefficient. As it can be seen, the calculated 
results exhibit some oscillations near the expansion head. 
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4. Inviscid Oblique Shock Reflection on a Plane Wall 


Numerical solutions of a regular shock reflection with an incident shock angle 
of 29 degrees and free stream Mach number of 2.9 have been obtained by solving the 
2D Euler equations. The calculations are started by assuming uniform flow with 
Mo„=2.9 everywhere except the top boundary, where the variables are consistently 
over specified from the jump conditions. The computational domain extends from 
x=0 to x=4 and from y=0 to y=l. Uniform grids are used to discretize the domain. 
Both coarse grid and fine grid calculations have been performed for the same CFL 
number. Table 2 lists the cases tested in this study. First of all, p > 0.5 seems to be a 
necessary condition for obtaining a converged solution. Secondly, the observation 
that, when (3=1.0, the coarse grid calculation without the use of flux limiter 
manages to reach a converged albeit smeared solution while its counter part fine 
grid calculation can not yield a converged solution, simply indicates that sharpened 
shock profile requires the use of flux limiter to maintain its accuracy and numerical 
stability. Thus, it is concluded that robust calculation of shocked flow requires the 
use of flux limiter in conjunction with a value of p > 0.5 . The pressure contours 
obtained from 122x42 grid points and with (3=1.0 are shown in Fig. 6(a). The 
comparison of computed and exact solutions along y=0.4878 is illustrated in Fig. 
6(b) and 6(c). Good results along other y=constant lines are also obtained. 

5. Supersonic Boundary Layer Flow Over a Flat Plate 

A Mach 2.2 laminar flow over an adiabatic flat plate is calculated by solving 
the 2D Navier-Stokes equations. The Reynolds number based on the free stream 
condition and a reference length of 0.08m is 9.8645xl0 4 . In the streamwise 
direction (- 0.19 < x < 2.00 with x=0 being the leading edge of the plate), 74 grid points 
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with a constant Ax=0.03 are used. In the direction normal to the wall (0 < y < 1.19212), 
62 grid points with the following distribution are used. 

Ayj = 1.5625xl0 4 , l<j<4 , 

Ayj = 1.1868A yj .i , 5 < j < 33 

Ayj = 3.75xl0‘ 2 , j < 33 (9) 

Several calculations have been performed to investigate the effects of (3 and 
the use of flux-limiter on the calculated boundary layer properties. All the 
calculations are carried out with a constant dimensionless time step At=0.01. Table 
3 shows the parameters of these cases and a short description of the most important 
observation. These results suggest that the use of flux-limiter does not degrade the 
computed boundary layer flow. Furthermore, the present Navier-Stokes solutions 
compare quite well with the corresponding solutions obtained by executing an 
existing 2D compressible boundary layer code [6], as illustrated in Fig. 7. These 
results are obtained with (3=1.0 and limiter on. For the evaluation of boundary 
layer properties such as the displacement and momentum thicknesses, the 
boundary layer edge is assumed to be at a point where the u-velocity is 99.5% of its 
free stream values. The profiles of u-velocity and temperature are normalized with 
free stream velocity and temperature respectively. They are then plotted against 
y/0 where 0 denotes the local momentum thickness. The skin friction coefficient is 
evaluated by using the first order backward differencing, and then plotted against 
Re x , which is the Reynolds number based on the local distance from the leading 
edge of the plate. 
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6. Interaction of Oblique Shock Wave and Laminar Boundary Layer. 


The main features of the oblique shock wave and laminar boundary layer 
interaction are well understood. A brief description of these features is given in Fig. 
8. Since this configuration involves shocks and separated flow, and carries no 
uncertainty associated with turbulence modeling, it is chosen here as a test case. 
The numerical solutions at the asymptotic steady state are compared with the 
experimental data reported in Ref. [7]. In addition, time-accurate calculations are 
performed to illustrate the transient development of this interaction. 

The flow parameters are such that M- = 2.2, Re = p„iux S h / |i~ = 9.8645xl0 4 , 
where X s h=0.08m is the distance between the leading edge and shock impingement 
point in the experiment. In addition, the shock incident angle is 30.027 degrees. 
Based on the studies discussed in Sections 4 and 5, the third-order-biased upwind 
(i.e., P=l) TVD scheme is used for the present calculations. The total number of grid 
points is 74x62, and they are distributed in two different ways, as depicted in Fig. 9. 
For a given grid distribution, a steady-state boundary layer flow is first established. 
Then, an oblique shock wave is imposed at the upper left corner of the inflow 
boundary and along the top boundary, where the variables are consistently over 
specified from the jump conditions. The subsequent transient development of the 
interacting process is followed by employing At=0. 01, until an asymptotic steady- 
state is reached. 

The grid distribution denoted by 'initial grids' is deduced from the triple-deck 
theory as described by Eq. (9) together with a constant Ax=0.03. The grid 
distribution denoted by 'adapted grids' is constructed by using the grid adaptation 
package TURBO AD [8]. More specifically, a steady state solution of the interaction 
is first obtained by using the initial grids. Using this grid distribution as a baseline 
and taking into account the large gradient regions in the flow solution, TURBO AD 
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enables a redistribution of grid points to improve the resolution in the shock and 
recirculation regions. Such an improvement is demonstrated in Fig. 10(a), which 
compares the measured wall pressure data and several computed results obtained 
from different grid distributions. It is noted here that, in Fig. 10(a), the pressure is 
normalized by Po, i.e., the minimum pressure just upstream of the interaction. The 
calculated skin friction distribution is shown in Fig. 10(b). Separation and 
reattachment points can be determined from the positions where the skin friction 
vanishes. Table 4 gives the locations of these points as determined from 
computations and experiment. The values determined from the adapted grids are 
in close agreement with the experimental values reported in Ref. [7]. Fig. 11 shows 
the measured and computed streamwise velocity profiles in the separated region. 
Since the results obtained from the adapted grids require further two-dimensional 
interpolations between grid points, the profiles obtained from the initial grid 
distribution are used for convenience. It is also noted here that, as discussed in Ref. 
[7], there are errors in the experimental velocity measurements due to the behavior 
of the seeding particles. 

The transient development of the interaction process is illustrated by the 
computed flow fields at three time stations, i.e., at the moment of shock 
impingement (t=3.3xl0' 4 seconds), the subsequent interaction in the upstream and 
downstream region (t=5.4xl0' 4 seconds), and the final asymptotic steady state. 
Figures 12 and 14(a), (b) describe the development of shock reflection and expansion 
fan in terms of the static pressure contours and wall pressure distribution 
(normalized by P„). Figures 13 and 14(c), (d) illustrate the emergence of a 
separation region and the thickening of the boundary layer in terms of the Mach 
number contours and skin friction distributions. It is noted here that the wall 
pressure and skin friction distributions are plotted in the computational space, i.e., 
the abscissa I is the streamwise grid index, and these results are obtained from the 
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adapted grid distribution. These figures, together with many others at different 
time stations, indicate that the entire process is essentially monotonic, i.e., no 
appreciable oscillation in flow patterns has been observed. 

7. Concluding Remarks 

Under the present effort, a TVD scheme has been developed and incorporated 
into an existing time-accurate high-resolution Navier-Stokes code. Test cases have 
been run to assess the overall accuracy and the robustness of the developed solution 
procedure. These numerical results compare well with corresponding analytic 
solution or experimental data. Its capability of accurately tracking transient 
movement of shock waves and contact discontinuities is demonstrated by 
calculating the shock tube problems. Its accuracy in calculating multi-dimensional 
shock waves and boundary layer flows is demonstrated by solving inviscid regular 
shock reflection and a supersonic boundary layer. It is then applied to compute the 
interaction of an oblique shock wave with a flat plate boundary layer, both the 
transient development and the asymptotic steady state of the interaction are 
presented to add some insight into this classical problem. 
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0 < x < 5 

5 < x < 10 
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u 

P 

P 
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P 

Sod's Problem 
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0 
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1 

0 
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Lax's Problem 

0.890 
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1 
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0 
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1 

0 

1 

Strong Shock 
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0 
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1 

0 
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Tablel. Normalized Initial Conditions for Shock Tube Cases 



3 

limiter 

Remark 

Coarse Grids (61x21) 

0.00 

on 

No solution 


0.10 

on 

No solution 


0.25 

on 



0.50 

on 



1.00 

on 



0.00 

off 

No solution 


0.25 

off 

No solution 


0.50 

off 



1.00 

off 


Fine Grids (122x42) 

0.25 

on 

No solution 


0.50 

on 



1.00 

on 



1.00 

off 

No solution 


Table 2. Cases Studied for Oblique Shock Wave Reflection in Inviscid Flow 
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3 

limiter 

Remark 

0.0 

off 

significant spatial oscillations in boundary layer 
thickness and skin friction 

0.5 

off 

no appreciable differences from P=1 and limiter-on 
case 

0.5 

on 

essentially the same as p=0.5 and limiter-off case 
except for some spatial oscillations in boundary layer 
thickness 

1.0 

off 

slightly higher wall temperature than the P=l, 
limiter-on case 

1.0 

on 

Results shown in Fig. 7 


Table 3. Cases Studied for Supersonic Laminar Boundary Layer Flow 



Separation (X/Xsh) 

Reattachment (X/Xsh) 

Present code (initial grids) 

0.745 

1.390 

Present code (adapted grids) 

0.795 

1.256 

Degrez' experiment 

0.78 

1.28 

Degrez' computation (M=2.15) 

0.79 

1.24 


Table 4. Separation and Reattachment Points in Shock Boundary Layer Interaction 
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Fig. 1. Sod's problem (J3=1.0) 
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Fig. 2. Lax's problem ((5=1.0) 
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Fig. 3. Large temperature ratio problem ((3=1.0) 
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Fig. 4. The appearance of expansion shock in computed flow (p=1.0 and a=0.0) 
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Fig. 5. Strong shock problem (p=1.0 and ct=1.5) 
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(a) pressure contours 



(b) pressure distribution (y=0.4878) 
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Fig. 6. Oblique shock reflection in inviscid flow 

(Mo. = 2.9, incoming shock angle=29°, P=1.0) 
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Fig. 7. Solutions for the supersonic laminar boundary layer over a flat plate 

(P=1.0 and limiter on) (a) normalized u-velocity profiles (b) normalized 
temperature profiles (c) skin friction coefficient (d) ratio of displacement 
thickness to momentum thickness 


19 


Page intentionally left blank 




21 


Figure 9. Computational domains 
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Fig. 10. Shock boundary layer interaction 
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Fig. 11. Streamwise velocity profiles (shock boundary layer interaction) 
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Figure 12. Pressure contours and wall pressure distributions at different times 
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Figure 13. Mach number contours ( Note that the x coordinate has been compressed ) 
and skin friction coefficient distributions at different times 
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Figure 14. Asymptotically converged solution for shock boundary layer interaction 
(Note that the x coordinate in Mach number contour has been compressed.) 
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